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ABSTRACT 

A shell/3D modeling technique was developed for which a local three-dimensional solid 
finite element model is used only in the immediate vicinity of the delamination front. The goal was 
to combine the accuracy of the full three-dimensional solution with the computational efficiency of 
a plate or shell finite element model. Multi-point constraints provided a kinematically compatible 
interface between the local three-dimensional model and the global structural model which has been 
meshed with plate or shell finite elements. Double Cantilever Beam (DCB), End Notched Flexure 
(ENF), and Single Leg Bending (SLB) specimens were analyzed first using three-dimensional 
finite element models to obtain reference solutions. Mixed mode strain energy release rate 
distributions were computed across the width of the specimens using the virtual crack closure 
technique. The analyses were repeated using the shell/3D technique to study the feasibility for pure 
mode I (DCB), mode II (ENF) and mixed mode I/II (SLB) cases. Specimens with a unidirectional 
layup and with a multidirectional layup where the delamination is located between two non-zero 
degree plies were simulated. For a local three-dimensional model, extending to a minimum of 
about three specimen thicknesses on either side of the delamination front, the results were in good 
agreement with mixed mode strain energy release rates obtained from computations where the 
entire specimen had been modeled with solid elements. For large built-up composite structures the 
shell/3D modeling technique offers a great potential for reducing the model size, since only a 
relatively small section in the vicinity of the delamination front needs to be modeled with solid 
elements. 
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INTRODUCTION 

One of the most common failure modes for composite structures is delamination. The 
remote loadings applied to composite components are typically resolved into interlaminar tension 
and shear stresses at discontinuities that create mixed-mode I and II delaminations. To characterize 
the onset and growth of these delaminations the use of fracture mechanics has become common 
practice over the past two decades [1-3]. The total strain energy release rate, G, the mode I 
component due to interlaminar tension, G,, the mode II component due to interlaminar sliding 
shear, G„, and the mode III component, G m , due to interlaminar scissoring shear, are calculated 
from continuum (2D) and solid (3D) finite element analyses using the virtual crack closure 
technique [4-8]. In order to predict delamination onset or growth, these calculated G components 
are compared to interlaminar fracture toughness properties measured over a range from pure mode 
I loading to pure mode II loading [9-12], 

Three-dimensional finite element models have been used to study the behavior of 
specimens used in fracture toughness testing [7,13-15], as well as the behavior of edge 
delaminations [1,16] and near-surface delaminations in composite laminates [17,18]. Since many 
layers of brick elements through the thickness are often necessary to model the individual plies, the 
size of finite element models required for accurate analyses may become prohibitively large. To 
improve computational efficiency, built-up structures are therefore traditionally modeled and 
analyzed using plate or shell finite elements. Computed mixed mode strain energy release rate 
components, however, depend on many variables such as element order and shear deformation 
assumptions, kinematic constraints in the neighborhood of the delamination front, and continuity 
of material properties and section stiffness in the vicinity of the debond when delaminations or 
debonds are modeled with plate or shell finite elements [7,19]. For example, in reference 19, mesh 
refinement studies showed that computed G h G„, and Gm did not converge when the structure 
above and below the plane of delamination was modeled with plate elements with different section 
properties (thickness or layup). A comparison of computed mixed mode strain energy release rates 
obtained from plate models with values computed from three-dimensional models showed 
differences in results near the free edges of the structure where the stress state is three-dimensional 
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[20]. These problems may be avoided by using three-dimensional models. Furthermore, three- 
dimensional analyses are required when matrix cracks and multiple delaminations need to be 
modeled at different ply interfaces. Three-dimensional analyses become necessary e.g. to analyze 
the skin/stringer debonding discussed in reference 21, where the failure at the flange tip is 
inherently three-dimensional as shown in Figure 1 . Matrix cracks and delaminations at different ply 
interfaces need to be modeled with solid elements. Therefore, methods based on three-dimensional 
modeling to calculate fracture parameters hi built-up structures need to be improved. 

The overall objective of the current work is to develop a shell/3D modeling technique for 
which a local solid finite element model is used only in the immediate vicinity of the delamination 
front and the remainder of the structure is modeled using plate or shell elements. The goal of the 
shell/3D technique is to combine the computational efficiency of a plate or shell finite element 
model with the accuracy of the full three-dimensional solution in the areas of interest. Multi-point 
constraints provide a kinematically compatible interface between the local three-dimensional model 
and the surrounding global structural model, which can be meshed with plate or shell finite 
elements. For large composite structures, the shell/3D modeling technique offers great potential for 
saving modeling and computational effort because only a relatively small section in the vicinity of 
the delamination front needs to modeled with solid elements as shown in Figure 2. A significant 
reduction hi model size can be expected compared to a full three-dimensional model. 

In the current investigation, the feasibility of the shell/3D technique proposed is studied for 
the pure mode I case, mode II case and a mixed mode I/II case. This is accomplished by using 
simple specimens like the double cantilever beam (DCB), end notched flexure (ENF), and single 
leg bending (SLB) specimens. First, all three specimens were modeled entirely with solid elements 
to validate the three-dimensional model, select suitable element types and an appropriate mesh size 
around the delamination front. These results were used as reference solutions for comparison with 
values obtained from the shell/3D technique. For each specimen, mixed mode strain energy release 
rate distributions were computed across the width from nonlinear finite element analyses using the 
virtual crack closure technique [3-5], The length of the local three-dimensional model around the 
delamination front was increased until the results computed were within 1% of mixed mode strain 
energy release rates obtained from computations where the entire specimen had been modeled with 
solid elements. 
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SPECIMEN DESCRIPTION 


For this investigation the double cantilever beam (DCB), the end notched flexure (ENF), 
and the single leg bending (SLB) specimens, as shown in Figure 3, were chosen to study the 
feasibility of the shell/3D technique for the pure mode I case, mode II case and a mixed mode I/II 
case, respectively, hi general DCB, ENF and mixed mode tests are performed on unidirectionally 
reinforced laminates, which means that delamination growth occurs at a [0/0] interface and crack 
propagation is parallel to the fibers. Although this unidirectional layup is desired for standard test 
methods to characterize fracture toughness, this kind of delamination growth will rarely occur in 
real structures. Previously, a number of combined experimental and numerical studies on 
unidirectional and multidirectional laminates have been performed where the critical strain energy 
release rates of various interfaces were evaluated under mode I, mode II and mixed-mode 
conditions [13, 14, 15, 22, 23]. Three different laminates were selected from these previous 
studies. The unidirectional layup [0] 32 was designated UD32, the unidirectional layup [0] 24 was 
designated UD24 and the multidirectional layup [ ± 3 ()/( )/- 3 ( )/( )/ 3 ( )/( ) 4 /3 0/ ()/- 3 ( )/0/- 3 ( )/ 3 0/ T - 
30/30/0/30/0/-30/0 4 /30/0/30/0/±30] was designated D±30. The arrow denotes the location of the 
delamination, which for all three laminates was located in the midplane. For interfacial 
delaminations between two orthotropic solids care must be exercised in interpreting the computed 
mixed inode energy release rates obtained from the virtual crack closure technique. This will be 
discussed in detail in the section on the SLB specimen with D±3() layup. The UD32 and D±30 
layup were made of C12K/R6376 graphite/epoxy and the UD24 layup was made of T300/1076 
graphite/epoxy. The material properties are given in Table 1 and the layup is summarized in 
Table 2. 


ANALYSIS FORMULATION 

FINITE ELEMENT ANALYSIS 

The goal of this investigation was to study the accuracy of the shell/3D modeling technique 
by comparing strain energy release rates computed using the shell/3D modeling technique to results 
obtained from full three-dimensional models. Therefore, all three specimens were first modeled 
entirely with solid elements. A typical three-dimensional finite element model of a specimen is 
shown in Figure 4(a). For the entire investigation, the ABAQUS® geometric nonlinear analysis 
procedure was used. To study the influence of element selection on the global load/deflection 
behavior and the computed mixed mode strain energy release rates, several three-dimensional solid 


- 4 - 



element types were used to model the specimens. The use of standard solid eight-noded brick 
element C3D8, incompatible mode element C3D8I, reduced integration element C3D8R as well as 
solid twenty-noded hexahedral elements C3D20 and reduced integration element C3D20R was 
studied. Shear locking is common in first-order, fully integrated elements, such as C3D8, that are 
subject to bending. The numerical formulation of tins element gives rise to shear strains that do not 
really exist. Therefore, these elements are too stiff in bending, and many elements over the 
thickness are required to obtain acceptable results. Elements where a lower-order, reduced 
integration is used to form the element stiffness such as the C3D8R and C3D20R elements usually 
provide more accurate results in bending and reduce running time. Incompatible mode elements, 
such as C3D8I are also recommended for bending and contact problems. In these elements internal 
deformation modes are added to the standard displacement modes of the element in order to 
eliminate the parasitic shear stresses that occur in bending [24], Interpenetration of the delaminated 
faces was prevented by using multi point constraints or contact elements [24], Results will be 
discussed in detail in the following chapter. 

The specimens with unidirectional layup were modeled by six elements through the 
specimen thickness as shown in the detail of Figure 4(b). For the specimens with D±30 layup, two 
plies on each side of the delamination were modeled individually using one element for each ply as 
shown in Figure 4(a). The adjacent four plies were modeled by one element with material 
properties smeared using the rule of mixtures [25] . The adjacent element extended over the four 0° 
plies. The six outermost plies were modeled by one element with smeared material properties. The 
delamination was modeled as a discrete discontinuity in the center of the specimen, with separate, 
unconnected nodes (with identical coordinates) on the upper and lower surfaces of the delaminated 
section. Referring to Figure 4, the specimens were divided into a center section of width, f, and a 
refined edge section, e, to capture local edge effects and steep gradients. These sections appear as 
dark areas in the full view of the specimen. Along the length of the model, a refined mesh of 
length, c, was used in the vicinity of the delamination front as shown in Figure 4(a). To study the 
influence of the mesh size around the delamination front on computed mixed mode strain energy 
release rates, the length, c, of the ref ined zone and the number of elements in the zone were varied, 
as will be discussed in detail later. 

A shell/3D model of a typical specimen is shown in Figure 4(b). The global section was 
modeled with ABAQUS® four-noded quadrilateral S4 type shell elements. The local three- 
dimensional section was modeled with ABAQUS® solid eight-noded C3D8I type elements. A 
combination of reduced integrated eight-noded quadrilateral shell elements S8R with solid twenty- 
noded hexahedral elements C3D20R was also studied. The transition from the global shell element 
model to the local three-dimensional model in the vicinity of the delamination front was 
accomplished by using multi-point constraint options given by ABAQUS® to enforce appropriate 
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translations and rotations at the shell-solid interface [24], The theory used for the multi-point 
constraint option assumes that the interface between the shell elements and solid elements is a 
surface containing the normals to the shell along the line of intersection of the meshes, so that the 
lines of nodes on the solid mesh side of the interface in the normal direction to the surface are 
straight lines. The nodes on the solid mesh side have the possibility of moving along the line and 
the line is allowed to change length, which means that there are no constraints in thickness 
direction [24], An improved coupling of the shell element model to the local three-dimensional 
model may be obtained by the use of special transition elements using formulations for the shell/3D 
transition based on a higher-order shell theory [26]. 

Along the length of the model a refined mesh of length c=5 mm with 10 elements was used 
for the UD24 layup, c= 3.0 mm, 12 elements for the D±30 layup, and c= 6.0 mm, 24 elements was 
used for the UD32 layup in the vicinity of the delamination front. These section lengths, c, had 
been selected in previous studies [14, 15, 23] and remained unchanged during the current 
investigations. To study the influence of the size of the local zone on computed mixed mode strain 
energy release rates, the total length of the local zone modeled with solid elements, d, was varied 
between 1 0 to 30 mm, in 5mm increments. 


VIRTUAL CRACK CLOSURE TECHNIQUE 


The Virtual Crack Closure Technique (VCCT) was used to calculate strain energy release 
rates [3-5], The mode I, mode II and mode III components of the strain energy release rate, G„ 
G,i and G,,,, were calculated for eight noded solid elements as shown in Figure 5 


G| = 
G„ = 

Gin = 


2AA 1 L 
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2AA 




— — Yu (v; 
2AA 1 L 


- W L/ ) 
- U L /) 
- V L /) 


with AA=Aab [5]. Here AA is the area virtually closed, Aa is the length of the elements at the 
delamination front and b is the width of the elements. For better identification in this and the 
following figures, columns are identified by capital letters and rows by small letters. Hence, X’ Li , 
Y’ u and Z’ u denote the forces at the delamination front in column L, row i, and u’ K „ v’ K/ and w’ KJ , are 

the relative displacements at the corresponding node row l behind the delamination front as shown 
in Figure 5. For geometrically nonlinear analysis, both forces and displacements were transformed 
into a local coordinate system (x’, y’), that defined the normal and tangential coordinate directions 
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at the delamination front in the defonned configuration. 

For twenty noded solid elements the equations to calculate the strain energy release rate 
components at the element comer nodes (location Li) as shown in Figure 6 are as follows [17] 



Here X’ Ki , Y’ Ki and Z’ Ki denote the forces at the delamination front in column K, row i, and u’ K( , v’ K/ 
and w’ w are the relative displacements at the corresponding column K, node row i behind the 
delamination front as shown in Figure 6. Similar definitions are applicable in column M for the 
forces at node row i and displacements at node row i and in column L for the forces at node row i 
and j and displacements at node row £ and m respectively. As mentioned previously, for 
geometrically nonlinear analysis, both forces and displacements were transformed into a local 
coordinate system (x’, y’), that defined the normal and tangential coordinate directions at the 
delamination front in the deformed configuration. The equations to calculate the strain energy 
release rate components at the mid side node (location Mi) as shown hi Figure 7 are as follows [17] 



Additional information with respect to the application of the VCCT and improved equations for 
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twenty noded solids are given in the literature [4-8, 17]. 

The total strain energy release rate, G T , was obtained by summing the individual mode 
components as 

Gr = G| +G| +Gm . 

The data required to perform the Virtual Crack Closure Technique were accessed from the 
ABAQUS® result file. The calculations were performed in a separate post-processing step using 
nodal displacements and nodal forces obtained from elements at the delamination front. 


ANALYSIS OF SPECIMENS WITH UNIDIRECTIONAL LAYUP 

Numerical validation of the finite element models and the post processing to compute the 
mixed mode strain energy release rates was performed using three-dimensional models of 
unidirectionally laminated DCB, ENF and SLB specimens. For each specimen, mixed mode strain 
energy release rate distributions were computed across the width from nonlinear finite element 
analyses using the virtual crack closure technique. Results were used as reference solutions for 
comparison with values obtained from the shell/3D technique. The analyses were then repeated 
using the shell/3D modeling technique to study the feasibility for pure mode I (DCB), mode II 
(ENF) and mixed mode I/II cases (SLB). 


COMPUTATION OF STRAIN ENERGY RELEASE RATES ACROSS A STRAIGHT 
DELAMINATION FRONT IN A DCB SPECIMEN WITH UNIDIRECTIONAL LAYUP 

For this investigation, the symmetry of the DCB specimen was taken into account and only 
one half of the specimen width B/2 was modeled as shown in Figure 8(a). The influence of 
element selection on the global load/deflection behavior and the computed mode I strain energy 
release rates was studied using ABAQUS® C3D8, C3D8I, C3D8R, C3D20 and C3D20R type 
elements to model the specimen. Along the length of the model, a refined mesh of length, c, was 
used in the vicinity of the delamination front. This section length, c, was varied to determine the 
length required for accurately computing the mixed mode strain energy release rates. For this 
study, the specimen was modeled with C3D8I elements and the number of elements, n, was varied 
accordingly to keep the element size constant as shown in Figures 8(a) to (c). Additionally, the 
influence of mesh size on computed mixed mode strain energy release rates was studied by keeping 
the length of the refined zone, c, constant, and increasing the number of elements, n, in this zone 
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as shown in Figures 9(a) to (c). The specimen was divided into a center section of width, f, and a 
refined edge section of width, e, to capture local edge effects and steep gradients. Section widths 
and mesh sizes used in the current investigation were taken from a previous study [23] where the 
effect of mesh size across the width on computed mixed mode energy release had been 
investigated. 

To study the influence of element selection on the global load/deflection behavior and the 
computed mixed mode strain energy release rates the mesh shown in Figure 8(b) was used for all 
element types. The load/deflection behavior of the specimen is shown in Figure 10 where the tip 
opening displacement 5, as shown in Figure 3(a), is plotted versus the applied load P. The 
deformation behavior computed using C3D8I, C3D20 and C3D20R elements matches previous 
results using solid twenty-noded hexahedral elements [23], the model therefore accurately captured 
the global response of the specimen. The classical eight noded brick element, C3D8, however, 
shows the tendency to lock, which means that an unnaturally stiff behavior of the structure is 
observed during computation. The eight noded brick element with reduced integration, C3D8R 
appears to model an excessively compliant structural behavior, hi Figure 1 1 the computed mode I 
strain energy release rate normalized with respect to the value from beam theory, 


G 


I, beam 



12- a 2 P 2 
B 2 h 3 ■ E-, 


is plotted versus the normalized width, y/B, of the specimen. Here, a denotes the delamination 
length, P the external load, B the specimen width, h the thickness of the cantilever arms as shown 
in Figure 3 and the modulus of elasticity. Strain energy release rates obtained from models 
using C3D8I, C3D20 and C3D20R elements are in excellent agreement with the values from the 
analyses using twenty noded solid elements [13]. The mode I strain energy release rate is fairly 
constant in the center part of the specimen progressively dropping towards the edges causing the 
straight front to grow into a curved front as explained in detail in the literature [8, 13, 27, 28]. As 
expected, the mode II and mode III strain energy release rates are computed to be nearly zero and 
hence are not shown. The results indicate that the post processing module used to compute the 
strain energy release rate using VCCT operates accurately. The model made of C3D8 elements 
yields results which do not correctly capture the drop of the mode I strain energy release rate 
towards the edge. Studies, where the model made of C3D8 elements was repeatedly refined in all 
three spatial directions verify that the load/deflection behavior and the strain energy release rates 
converge to the results reported for the higher order elements. The mixed mode strain energy 
release rates computed from the model made of C3D8R elements are noticeably higher across the 
entire width compared to the other results. The response of specimens modeled with C3D8R 
elements was not investigated any further. Compared to the model made of C3D20 elements the 
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models made of C3D8I and C3D20R yield nearly the same results and require less computation 
time. Therefore these element types were chosen for the following studies. 

The influence of the length, c, of the refined section around the delamination tip (Figure 8) 
on the computed mode I strain energy release rate distribution across the width of a DCB specimen 
was studied next for a model made of C3D8I type elements. As the length of the refined section 
was increased, the element length was kept constant at Aa=0.5 mm. The plot in Figure 12 indicates 
that the length, c, has only a small influence on the computed strain energy release rate distribution. 
In order to be consistent with previous studies [13, 23], a length of c= 5 mm was chosen for the 
local section of the shell/3D model. Additionally, the influence of mesh size was studied. The mesh 
size is defined as the length of the elements in the refined section, c, which is identical to the 
length, Aa, of the elements at the delamination tip. The influence on the mode I strain energy 
release rate distribution across the width is moderate as shown in Figure 13 and only very long 
elements (n=2, Aa =5 mm) need to be avoided. Hence, ten elements (n=10) were chosen for the 
local three-dimensional section of the shell/3D model, which leads to an element length of Aa=0.5 
mm at the delamination tip. 

A shell/3D model of a DCB specimen is shown in Figure 14(a). The global section was 
modeled with S8R type shell elements. The local three-dimensional section was modeled with solid 
C3D20R type elements. A combination of quadrilateral shell elements S4 with solid eight-noded 
elements C3D8I, was also studied. Along the length a refined mesh of length c=5 mm with 10 
elements was used in the vicinity of the delamination front. This section length, c, was kept 
constant during the entire investigation. To study the influence of the length of the local zone on 
computed mixed mode strain energy release rates, the total length of the local zone modeled with 
solid elements, d, was varied (d=10, 15, 20, 25, 30 mm) as shown in Figures 14 (a) to (f). 

The computed strain energy release rate distributions across the width of the specimen are 
shown in Figure 15 for the combination of twenty noded brick elements in the local three- 
dimensional model with eight noded shell elements in the global model. As expected, the mode II 
and mode III strain energy release rates are computed to be nearly zero and hence are not shown. 
With increasing length of the local three-dimensional model, d, computed results from the shell/3D 
model converge to the solution obtained from a full three-dimensional model. For a local three- 
dimensional model extending to a minimum of about three specimen thicknesses in front and 
behind the delamination front (d/2h=6.67), the results were within 1% of the mode I strain energy 
release rates obtained from computations where the entire specimen had been modeled with solid 
elements. The shell/3D model is capable of accurately simulating the anticlastic bending effect that 
causes the strain energy release rate to be highest in the center of the specimen and lowest at its 
edges. Results for the combination of eight noded brick elements with four noded shell elements 
are shown in Figure 16. When compared to the model with twenty noded solids and eight noded 
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shell elements, the combination of eight noded brick elements with four noded shell elements 
yields identical results and provides a reduced model size. Therefore, the combination of eight 
noded brick elements with four noded shell elements was used for the remainder of this study. 


COMPUTATION OF STRAIN ENERGY RELEASE RATES ACROSS A STRAIGHT 
DELAMINATION FRONT IN AN ENF SPECIMEN WITH UNIDIRECTIONAL LAYUP 

For this investigation, the symmetry of the ENF specimen was taken into account and only 
one half of the specimen width B/2 was modeled as shown in Figure 17. The influence of element 
selection on the computed mode II strain energy release rates was studied using ABAQUS® C3D8, 
C3D8I, C3D8R, C3D20 and C3D20R type elements to model the specimen. Along the length of 
the model a section c=10 mm was used for the refined mesh in the vicinity of the delamination 
front. Twenty elements were used yielding an element length Aa=c/n=0.5 mm. This mesh size was 
found suitable during the initial investigation of the DCB specimen described in the previous 
section. Interpenetration of the cantilever arms was first prevented by using contact elements [24], 
A previous study [23] showed that penetration of the anus could be prevented by introducing multi 
point constraints in the plane of delamination only along a string of nodes above the left support as 
schematically shown hi Figure 18. The use of multi point constraints appears advantageous as less 
modeling effort is required and a computationally expensive contact analysis is avoided. Therefore, 
the influence of the multi point constraint technique on the computed mode II and mode III strain 
energy release rates was also studied. 

In Figures 19 aud 20 the computed mode II and mode III strain energy release rates 
nonnalized with the reference value, G llbeam , from classical beam theory (not accounting for 
transverse shear) 


l,bea 


d a ) 


9 a 2 P 2 
16 B 2 - h 3 • E-, ’ 


are plotted versus the nonnalized width, y/B, of the specimen. Computed mode II and mode III 
strain energy release rates obtained from models using C3D8I, C3D20 and C3D20R elements are 
in excellent agreement with the values from previous analyses using twenty noded solid elements 
[13]. The mode II strain energy release rate is fairly constant across almost the entire width of the 
specimen, peaking in the immediate vicinity of the edges. The mode III contribution is zero in the 
center of the specimen peaking to about only 5% of G M beam at the edges. The computed G values 
are nearly zero and therefore are not shown. The model made of C3D8 elements yields results 
which do not correctly capture the mode II and III distribution. The mode II values computed from 
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the model made of C3D8R elements appeal* excessively high across the entire width. These 
observations support the results obtained from the study of the UD24 DCB specimen discussed 
above. Models made of C3D8 and C3D8R elements were not investigated any further. Compared 
to the model made of C3D20 elements the models made of C3D8I elements yield nearly the same 
results and provide a reduced model size. Therefore this element type was chosen for the following 
studies. 

Computed mode II and mode III strahi energy release rates obtained from models where the 
penetration of the cantilever arms was prevented by multi point constraints or contact analysis are 
shown in Figure 21. The results shown are in almost exact agreement. Therefore, the technique 
using the multi point constraints to avoid penetration was used for the remainder of this study. 

A shell/3D model of an ENF specimen is shown in Figure 22. The global section was 
modeled with quadrilateral shell elements S4. The local three-dimensional section was modeled 
with eight-noded C3D8I type elements. Along the length, a refined mesh of length c=5 mm with 10 
elements was used in the vicinity of the delamination front. As discussed in the study of the DCB 
specimen above, the section length, c, was kept constant during the entire investigation. To study 
the influence of the length of the local zone on computed mixed mode strain energy release rates, 
the total length of the local zone modeled with solid elements, d, was varied (d=1 0, 1 5, 20, 25, 30 
mm) as shown in Figures 14 (a) to (f). 

The computed mode II and III strain energy release rate distributions across the width of 
the specimen are shown in Figures 23 and 24. As expected, the mode I strain energy release rate is 
computed to be nearly zero and hence is not shown. With increasing length of the local three- 
dimensional model, d, computed results from the shell/3D model converge to the solution obtained 
from a full three-dimensional model. For a local three-dimensional model extending to a minimum 
of about three specimen thicknesses in front and behind the delamination front (d/2h=6.67), the 
results were within 1% of the mode II strain energy release rates obtained from computations 
where the entire specimen had been modeled with solid elements. The results obtained from all 
shell/3D models are in excellent agreement with the mode III strain energy release rate obtained 
from computations where the entire specimen had been modeled with solid elements as shown in 
Figure 24. For a unidirectional layup, however, the mode III contribution is very small. 
Therefore, the influence of the size of the local three-dimensional model on mode II and mode III 
separation needs to verified for a case where the mode III contribution is more obvious and the 
mode II contribution is less dominant. This will be discussed in the section on the ENF specimen 
with multidirectional layup. 
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COMPUTATION OF STRAIN ENERGY RELEASE RATES ACROSS A STRAIGHT 
DELAMINATION FRONT IN A SLB SPECIMEN WITH UNIDIRECTIONAL LAYUP 

The single leg bending (SLB) specimen, as shown in Figure 3(c), was introduced for the 
determination of fracture toughness as a function of mixed mode I/II ratio [15]. This test may be 
performed in a standard three point bending fixture such as that used for the ENF test. By varying 
the relative thickness of the delaminated regions (t, and tg) various mode mixities may be achieved. 
The test is of particular interest because compliance calibration can be used to accurately determine 
the critical strain energy release rate [15]. This type of specimen was selected for this study to 
verify that the shell/3D modeling technique is also capable of accurately simulating the mixed mode 
I/II case. Mixed mode strain energy release rates which served as reference solutions had been 
computed in a previous study using three-dimensional FE models [23] . 

For this investigation, the symmetry of the SLB specimen was taken into account and only 
one half of the specimen width B/2 was modeled as shown in Figure 25. The finite element model 
is basically identical to the one used for the ENF and DCB specimens discussed in the previous 
sections, except the boundary conditions were modified by omitting the suppressions simulating 
the lower support pin of the ENF test. The influence of element selection on the computed mixed 
mode strain energy release rates was studied usmg ABAQUS® C3D8, C3D8I, C3D8R, C3D20 and 
C3D20R type elements to model the specimen. Along the length of the model, a section c=6 mm 
was used for the refined mesh in the vicinity of the delamination front. Twenty four elements were 
used yielding an element length Aa=c/n=0.25 mm. This mesh size was found suitable during the 
initial investigation of the DCB specimen described in the previous section. 

In Figures 26 to 28 the computed mode I, II and mode III strain energy release rates are 
plotted versus the normalized width, y/B, of the specimen. Computed strain energy release rates 
obtained from models using C3D8I, C3D20 and C3D20R elements are in excellent agreement with 
the values from a previous analysis using continuum based shell elements [23, 29] . As shown in 
Figure 26, the mode I strain energy release rate is fairly constant in the center part of the specimen 
progressively dropping towards the edges as previously discussed for the DCB specimen. The 
model made of C3D8 elements yields results which do not correctly capture the drop of the mode I 
strain energy release rate towards the edge. The values computed from the model made of C3D8R 
elements appear excessively high. The mode II strain energy release rate as shown in Figure 27 is 
fairly constant across almost the entire width of the specimen, peaking in the immediate vicinity of 
the edges, which was discussed earlier in the section about the ENF specimen. As shown in 
Figure 28, the mode III contribution is zero in the center of the specimen peaking to about only 
8% of G m at the edges. The model made of C3D8 elements yields results which do not correctly 
capture the mode II and III distribution. The mode II values computed from the model made of 
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C3D8R elements appeal - excessively high across the entire width of the specimen. These 
observations support the results obtained from the study of the UD24 DCB and ENF specimens 
discussed previously. Models made of C3D8 and C3D8R elements were not investigated any 
further. Compared to the model made of C3D20 elements the models made of C3D8I elements 
yield nearly the same results and provide a reduced model size. Therefore this element type was 
chosen for the following studies. 

A shell/3D model of a SLB specimen is shown in Figure 29. The global section was 
modeled with quadrilateral shell elements S4. The local three-dimensional section was modeled 
with eight-noded C3D8I type elements. Along the length a refined mesh of length c=6 mm with 24 
elements was used in the vicinity of the delamination front. As discussed in the study of the DCB 
and ENF specimens above, the section length, c, was kept constant during the entire investigation. 
To study the influence of the length of the local zone on computed mixed mode strain energy 
release rates, the total length of the local zone modeled with solid elements, d, was varied (d=10, 
15, 20, 25, 30 mm) as shown hi Figures 14 (a) to (f). 

The computed mode I, II and III strain energy release rate distributions across the width of 
the specimen are shown in Figures 30 to 32. With increasing length of the local three-dimensional 
model, d, computed results from the shell/3D model converge to the solution obtained from a full 
three-dimensional model. For a local three-dimensional model extending to a minimum of about 
three specimen thicknesses in front and behind the delamination front (d/2h=6. 1 6), the results were 
within 1% of the mode I and II strain energy release rates obtained from computations where the 
entire specimen had been modeled with solid elements. As shown in Figure 32, the results 
obtained from all shell/3 D models are in excellent agreement with the mode III strain energy release 
rate obtained from computations where the entire specimen had been modeled with solid elements. 
For unidirectional layup, however, the mode III contribution is very small. Therefore, the 
influence of the size of the local three-dimensional model on in plane shear mode II and mode III 
separation needs to verified for a case where the mode III contribution is more obvious and the 
mode II contribution is less dominant. This will be discussed in the sections on the ENF and SFB 
specimens with multidirectional layup. 
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ANALYSIS OF SPECIMENS WITH A MULTIDIRECTIONAL LAYUP 


In general DCB, ENF and mixed mode tests are performed on unidirectionally reinforced 
laminates, which means that delamination growth occurs at a [0/0] interface and crack propagation 
is parallel to the fibers. Although this unidirectional layup is desired for standard test methods to 
characterize fracture toughness, this kind of delamination growth will rarely occur in real 
structures. Previously, a number of combined experimental and numerical studies on specimens 
with multidirectional layup have been performed where the critical strain energy release rates of 
various interfaces were evaluated under inode I, mode II and mixed-mode conditions [14, 15, 23]. 
In this study, DCB, ENF and SLB specimens with a multidirectional layup were first modeled 
entirely with solid elements to validate the three-dimensional model and select an appropriate mesh 
size around the delamination front. For each specimen type, mixed mode strain energy release rate 
distributions were computed across the width from nonlinear finite element analyses using the 
virtual crack closure technique. Results were used as reference solutions for comparison with 
values obtained from the shell/3D technique. The analyses were then repeated using the shell/3D 
modeling technique to study the feasibility for pure mode I (DCB), mode II (ENF) and mixed 
mode I/II cases (SLB). 


COMPUTATION OF STRAIN ENERGY RELEASE RATES ACROSS A STRAIGHT 
DELAMINATION FRONT IN A SLB-TYPE SPECIMEN 

Previous investigations have shown that care must be exercised in interpreting the values 
for G„ G,i and G m obtained using the virtual crack closure technique for interfacial delaminations 
between two orthotropic solids [30, 31], Mathematical solutions of the near crack tip field indicate 
that stresses start to oscillate in the immediate vicinity of the tip when crack growth occurs at 
interfaces between materials with dissimilar properties, hi the current investigation, this 
phenomenon has to be considered as the delamination growth occurs at a +30%30° interface. 
Therefore, the mixed mode SLB specimen was studied first and an appropriate mesh size was 
determined which was then also used for the models of the D±30 DCB and ENF specimens. 

For the investigation of the D±30 SLB specimen, the finite element model shown in 
Figure 33 was used. The model was made of eight noded ABAQUS® C3D8I elements. It had been 
shown above that compared to the model made of C3D20 elements the models made of C3D8I 
yield nearly the same results and provide a reduced model size. Therefore these element types were 
chosen for this and the following studies. The specimen was divided into a center section of width, 
f, and a refined edge section of width, e, to capture local edge effects and steep gradients. Section 
widths and mesh sizes used in the current investigation were taken from a previous study [23] 
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where the effect of mesh size across the width on computed mixed mode energy release had been 
investigated. Along the length of the model a refined mesh of length c was used in the vicinity of 
the delamination front. This section length, c=3.0 mm, was kept constant. The number of elements, 
n, was varied to study the influence of mesh size on computed mixed mode strain energy release 
rates. 

First, the influence of mesh size was studied. The mesh size equals the length of the 
elements in the refined section, which is identical to the length of the elements at the delamination 
tip. As shown in Figure 34, the zone with a constant G, distribution in the center becomes smaller 
compared to the UD32 case and the drop towards the edges is more pronounced. The drop is 
caused by increased anticlastic bending due to the lower values of bending rigidities in the 
individual arms. The influence of mesh refinement on the mode I strain energy release rate 
distribution across the width is moderate and only very long elements (n=3, Aa=c/n= 1 mm) need to 
be avoided. This is confirmed by the mode II and mode III distributions as shown in Figures 35 
and 36 where the mode II strain energy release rate is fairly constant across almost the entire width 
of the specimen and peaks near the edges accompanied by local mode III contribution. Compared 
to the UD32 layup these peaks become more visible for specimens with the D±30 layup caused by 
increased anticlastic bending. The distribution of the mixed mode ratio G,/G M is shown in Figure 
37. For the range studied (n=3 up to 48), there is only a small dependence of computed mixed 
mode ratio on element size Aa=c/n. Hence, twelve elements were chosen (n=1 2) in the refined 
section (c=3 mm). For the delamination in the +307-30° ply interface, the element length was 
therefore chosen to be Aa=c/n=0.25 mm and this element length was used consistently during the 
entire investigation. 

Second, the computed mode I, II and mode III strain energy release rates as shown in 
Figures 38 to 40 were compared with values from previous analyses using layered, continuum 
based shell elements [23, 29], The good agreement of the results indicates that the model was set 
up appropriately and the post processing module to compute the strain energy release rate using 
VCCT operates accurately. For comparison, mixed mode strain energy release rates were 
computed from models where the local penetration of the cantilever anus at the specimen edge near 
the delamination front was prevented by contact analysis. The results included in Figure 38 to 40 
are in almost exact agreement with results from simple analyses, where the penetration was not 
prevented. It was therefore chosen not to prevent the penetration and thus avoid the complicated 
contact analysis for the remainder of this study. 

A shell/3D model of a SLB D±30 specimen is shown in Figure 41. The global section was 
modeled with quadrilateral S4 type shell elements. The local three-dimensional section was 
modeled with eight-noded C3D8I type elements. Along the length a refined mesh of length c=3 mm 
with 12 elements was used in the vicinity of the delamination front. As discussed earlier, the 
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section length, c, was kept constant during the entire investigation. To study the influence of the 
length of the local zone on computed mixed mode strain energy release rates, the total length of the 
local zone modeled with solid elements, d, was varied (d=1 0, 1 5, 20, 25, 30 mm) as shown in 
Figures 14 (a) to (f). 

The computed mode I, II and III strain energy release rate distributions across the width of 
the specimen are shown in Figures 42 to 44. With increasing length of the local three-dimensional 
model, d, computed results from the shell/3D model converge to the solution obtained from a full 
three-dimensional model. For a local three-dimensional model extending to a minimum of about 
three specimen thicknesses in front and behind the delamination front (d/2h=6. 1 6), the results were 
within 1% of the mode I and II strain energy release rates obtained from computations where the 
entire specimen had been modeled with solid elements. As shown in Figure 43, the computed 
mode II contribution is constant across the center of the specimen, with larger values near the 
edges compared to the results obtained for the UD32 laynp. For the D±30 layup also a 
considerable amount of mode III is present due to the increased anticlastic bending for this layup as 
shown in Figure 44. As before, the results from the shell/3D technique were within 1% of the 
reference solution obtained from computations where the entire specimen had been modeled with 
solid elements. These results indicate that the shell/3D technique is capable of accurately simulating 
the increased anticlastic bending effect due to the lower values of bending rigidities in the 
individual arms for this layup which causes the mode II and III strain energy release rate to be 
higher towards the free edges. 


COMPUTATION OF STRAIN ENERGY RELEASE RATES ACROSS A STRAIGHT 
DELAMINATION FRONT IN A DCB SPECIMEN 

For this investigation, the entire width B of the DCB specimen was modeled as shown in 
Figure 45. The eight noded brick element ABAQUS® C3D8I was used for the simulation. Along 
the length of the model, a refined mesh of length c= 3 mm was used in the vicinity of the 
delamination front. The refined section was divided into n= 1 2 number of elements, which was 
found to yield a reasonable mesh Aa=c/n=0.25 mm as discussed for the SLB specimen in the 
previous section. The specimen was divided into a center section of width, f, and a refined edge 
section of width, e, to capture local edge effects and steep gradients. Section widths and mesh 
sizes used in the current investigation were taken from a previous study [23] where the effect of 
mesh size across the width on computed mixed mode energy release had been investigated. 

The computed mode I strain energy release rate distribution as shown in Figure 46 was 
compared with values from a previous analysis using layered, continuum based shell elements [23, 
29]. The mode I strain energy release rate is fairly constant in the center part of the specimen 
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progressively dropping towards the edges causing the straight front to grow into a curved front. 
Basically the distributions are similar to those computed for the UD24-layup. For specimens with 
multidirectional layup the zone with a constant mode I distribution in the center becomes smaller 
and the drop towards the edges is more pronounced. This phenomenon is caused by the smaller 
bending rigidities in the individual arms of the specimens and has been the subject of detailed 
experimental and analytical investigations [14, 23]. As expected, the mode II and mode III strain 
energy release rates are computed to be nearly zero and hence are not shown. The good agreement 
of the results indicates that the model was set up appropriately and the post processing module to 
compute the strain energy release rate using VCCT operates accurately. The computed mode I 
strain energy release rate obtained from a model where the penetration of the cantilever arms was 
prevented near the delamination front using contact analysis was also included in Figure 46 for 
comparison. The distribution shown is in almost exact agreement with the distribution obtained 
from a simple analysis, where the penetration was not prevented. It was therefore chosen not to 
enforce contact and thus avoid the complicated contact analysis for the remainder of this study. 

A shell/3D model of a DCB specimen is shown in Figure 47. The global section was 
modeled with S4 type shell elements. The local three-dimensional section was modeled with 
C3D8I type solid elements. Along the length a refined mesh of length c=3 mm with 12 elements 
was used in the vicinity of the delamination front, which is identical to the refined mesh used for 
the full three-dimensional model discussed above. This section length, c, was kept constant during 
the entire investigation. To study the influence of the length of the local zone on computed mixed 
mode strain energy release rates, the total length of the local zone modeled with solid elements, d, 
was varied (d=10, 15, 20, 25, 30 mm) as shown in Figures 14 (a) to (f). 

The computed strain energy release rate distributions across the normalized width of the 
specimen are shown in Figure 48. As expected, the mode II and mode III strain energy release 
rates are computed to be nearly zero and hence are not shown. With increasing length of the local 
three-dimensional model, d, computed results from the shell/3D model converge to the solution 
obtained from a full three-dimensional model. For a local three-dimensional model extending to a 
minimum of about three specimen thicknesses in front and behind the delamination front 
(d/2h=6.16), the results were within 1% of the mode I strain energy release rates obtained from 
computations where the entire specimen had been modeled with solid elements. As shown in 
Figure 48, the zone with a constant G, distribution in the center becomes smaller compared to the 
UD24 case and the drop towards the edges is more pronounced. The drop is caused by increased 
anticlastic bending due to the lower values of bending rigidities in the individual arms for this 
layup. The good agreement with results obtained from full three-dimensional models suggests that 
these effects are accurately simulated by the shell/3D model. 
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COMPUTATION OF STRAIN ENERGY RELEASE RATES ACROSS A STRAIGHT 
DELAMINATION FRONT IN AN ENF SPECIMEN 

For the mode II ENF tests, references 14 and 23 show that the mode II strain energy 
release rate is fairly constant across almost the entire width of the specimen, peaking in the 
immediate vicinity of the edges and accompanied by local mode III contributions in the same areas. 
These peaks become more visible for specimens with multidirectional layup [14, 23]. Therefore an 
ENF specimen with D±30 layup was selected for this study to verify the accuracy of the shell/3D 
technique in simulating this local mixed mode case near the edge of the specimen. For this 
investigation, the entire width B of the ENF specimen was modeled as shown in Figure 49. Along 
the length of the model a refined mesh of length c= 3 mm was used in the vicinity of the 
delamination front. The refined section was divided into n= 1 2 number of elements, which was 
found to yield a reasonable mesh with element size Aa=c/n=0.25 mm as discussed for the SLB 
specimen in the previous section. The specimen was divided into a center section of width f and a 
refined edge section of width e to capture local edge effects and steep gradients. Section widths 
and mesh sizes used in the current investigation were taken from a previous study [23] where the 
effect of mesh size across the width on computed mixed mode energy release had been 
investigated. Interpenetration of the cantilever anus was first prevented by using contact elements 
[24] . Earlier studies [23] showed, that penetration of the arms could be prevented by introducing 
multi point constraints in the plane of delamination only along a string of nodes above the left 
support as schematically shown in Figure 18. The use of multi point constraints appears 
advantageous as less modeling effort is required and a computationally expensive contact analysis 
is avoided. Therefore, the influence of the multi point constraint technique on the computed mode 
II and mode III strain energy release rates was also studied. 

In Figures 50 and 5 1 the computed mode II and mode III strain energy release rates are 
plotted versus the normalized width, y/B, of the specimen. The results are in good agreement with 
the distribution from a previous analysis using layered, continuum based shell elements [23, 29]. 
The mode II strain energy release rate is fairly constant across almost the entire width of the 
specimen, peaking near the edges and accompanied by local mode III contributions in the same 
area. Compared to the UD32 layup these peaks become more visible for specimens with D±30 
layup caused by increased anticlastic bending effect due to the lower values of bending rigidities in 
the individual arms for this layup. The computed G, values are nearly zero and therefore are not 
shown. Computed mode II and mode III strain energy release rates obtained from models where 
the penetration of the cantilever arms was prevented by multi point constraints were included in 
Figures 50 and 5 1 . The results shown are in almost exact agreement with the values obtained from 
a contact analysis. Therefore, the technique using the multi point constraints to avoid penetration 
was used for the remainder of this study. 
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A shell/3D model of an ENF specimen is shown in Figure 52. The global section was 
modeled with quadrilateral S4 type shell elements. The local three-dimensional section was 
modeled with eight-noded C3D8I type elements. Along the length a refined mesh of length c=3 mm 
with 12 elements was used in the vicinity of the delamination front. As discussed in the study of 
the DCB specimen above, the section length, c, was kept constant during the entire investigation. 
To study the influence of the length of the local zone on computed mixed mode strain energy 
release rates, the total length of the local zone modeled with solid elements, d, was varied (d=10, 
15, 20, 25, 30 mm) as shown in Figures 14 (a) to (f). 

The computed mode II and III strain energy release rate distributions across the normalized 
width of the specimen are shown in Figures 53 and 54. As expected, the mode I strain energy 
release rates is computed to be nearly zero and hence is not shown. With increasing length of the 
local three-dimensional model, d, computed results from the shell/3D model converge to the 
solution obtained from a full three-dimensional model. For a local three-dimensional model 
extending to a minimum of about three specimen thicknesses in front and behind the delamination 
front (d/2h=6.16), the results were within 1% of the mode II and III strain energy release rates 
obtained from computations where the entire specimen had been modeled with solid elements. As 
shown in Figure 53, the computed inode II contribution is constant across the center of the 
specimen, with larger values near the edges compared to the results obtained for the UD24 layup. 
For the D±30 layup also a considerable amount of mode III is present due to the increased 
anticlastic bending for this layup as shown in Figure 53. These results indicate that the shell/3D 
technique is capable of accurately simulating the increased anticlastic bending effect due to the 
lower values of bending rigidities in the individual arms for this layup which causes the mode II 
and III strain energy release rate to be higher towards the free edges. 
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CONCLUDING REMARKS 


A shell/3D modeling technique was presented for the analysis of composite laminates with 
delaminations. The individual mode and total strain energy release rates along the delamination 
front were evaluated. In this analysis, a local solid finite element model was used only in the 
immediate vicinity of the delamination front. The goal was to combine the accuracy of the full 
three-dimensional solution with the computational efficiency of a plate or shell finite element 
model. Multi-point constraints provided a kinematically compatible interface between the local 
three-dimensional model and the global structural model which was meshed with shell finite 
elements. 

For DCB, ENF, and SLB specimens, mixed mode strain energy release rate distributions 
were computed across the width from nonlinear finite element analyses using the virtual crack 
closure technique. This served to study the feasibility of the proposed shell/3D modeling technique 
for the pure mode I case (DCB), inode II case (ENF) and a mixed mode I/II case (SLB). 
Specimens with a unidirectional layup, for which the delamination is located between two 0° plies, 
as well as a multidirectional layup were simulated. First, all three specimens were modeled entirely 
with solid elements to validate the three-dimensional model, to select suitable element types and 
appropriate mesh size around the delamination front, and to check the need for contact analysis to 
prevent the interpenetration of the delaminated surfaces. Results were used as reference solutions 
for comparison with values obtained from the shell/3D technique. The geometrically nonlinear 
solution option of the ABAQUS® finite element code was used for the entire investigation. For 
each specimen, mixed mode strain energy release rate distributions were computed across the 
width from nonlinear finite element analyses using the virtual crack closure technique. 

The computation of mixed mode strain energy release rates is most critical for interfacial 
delaminations between two different orthotropic solids. Therefore, general recommendations for 
the selection of element types and appropriate mesh size around the delamination front may be 
taken from the results obtained from the specimens where the delamination is located between two 
non-zero plies. Compared to earlier studies, the models made of solid twenty-noded hexahedral 
elements (ABAQUS® types C3D20 and C3D20R with reduced integration) and solid eight-noded 
incompatible mode elements (ABAQUS®, type C3D8I) yield excellent results. The mesh 
refinement study showed that only a section of about 1 mm on either side of the delamination front 
needs to be refined. The influence of mesh size on the computed mixed mode ratio was negligible 
for elements lengths between Aa=0.5 mm (four ply thicknesses) and Aa=0.0625 mm (half a ply 
thickness). For the ENF specimen, it was found that instead of a complicated contact analysis, 
penetration of the arms could be prevented by introducing multi-point constraints between the 
delaminated surfaces just above the left support without compromising the accuracy of the 
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computed results. For the DCB and SLB specimens the influence of interpenetration was 
negligible. 

For the current investigation, two shell/3D combinations were studied: eight noded solid 
elements in the local section combined with four noded shell elements in the global section of the 
model and twenty-noded solid elements in the local section combined with eight-noded shell 
elements in the global section. The shell elements were connected to the local three-dimensional 
model using multi-point constraints to enforce appropriate translations and rotations. An overview 
of all element types, shell/3D combinations and mesh refinements investigated is given in Table 3 
for all specimen types. 

Finite element analyses showed that the accuracy achieved depends on the size of the local 
area. With increasing size of the local three-dimensional model, the computed results converged 
towards the strain energy release rates obtained from full three-dimensional finite element analysis 
for all specimens, layups and element types simulated. The results were in good agreement with 
the reference solution once the local zone was extended to about three times the specimen thickness 
in front and behind the delamination front. 

For large composite structures the shell/3D modeling technique offers great potential for 
reducing the model size because only a relatively small section in the vicinity of the delamination 
front needs to modeled with solid elements. A significant reduction in model size can be expected 
compared to a full three-dimensional model. In the current investigation, the application of the 
shell/3D technique reduced the number of degrees of freedom by about 35% compared to a full 
three-dimensional model for all three specimen types. Existing plate models may be modified to 
shell/3D models, which is a considerable advantage compared to the creation of an entirely new 
three-dimensional finite element model. 
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TABLES 


TABLE 1. MATERIAL PROPERTIES. 


T300/1076 Unidirectional Graphite/Epoxy Prepreg [16] 

E n = 139.4 GPa 

E 2 2 = 10.16 GPa 

£ 33 = 10.16 GPa 

V 12 = 0.30 

v 13 =0.30 

V 23 = 0.436 

G\ 2 = 4.6 GPa 

G 13 = 4.6 GPa 

G 23 = 3.54 GPa 

C12K/R6376 Unidirectional Graphite/Epoxy Prepreg [16] 

En = 146.9 GPa 

E 2 2 = 10.6 GPa 

£33 = 10.6 GPa 

V 12 = 0.33 

V 13 = 0.33 

V 23 = 0.33 

Gi 2 = 5.45 GPa 

G 13 = 5.45 GPa 

G 23 = 3.99 GPa 


TABLE 2. STACKING SEQUENCE. 

Layup-ID 

Stacking Sequence 

Material 

UD24 

[ 0] 24 

T300/1076 

UD32 

[ 0] 32 

C12K/R6376 

D±30 

[ ± 3 0/0/-3 0/0/3 0/0 4/3 0/0/- 3 0/0/- 3 0/3 0/ 1 - 3 0/3 0/ 3 0/0/30/0/ - 3 0/0 4/- 3 0/0/30/0/± 3 0] 

C12K/R6376 
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TABLE 3. OVERVIEW OF FINITE ELEMENT ANALYSES. 



DCB 

ENF 

SLB 

DCB 

ENF 

SLB 

Remarks 


UD24 

UD24 

UD32 

D±30 

D±30 

D±30 


element 







C3D8I, C3D20 and C3D20R 

type 1 







elements yielded almost identical 

C3D8 

+ 

+ 

+ 




results. 

C3D8I 

+ 

+ 

+ 




C3D8I was found to be 

C3D8R 

+ 

+ 

+ 




computationally most efficient 

C3D20 

+ 

+ 

+ 





C3D20R 

+ 

+ 

+ 





element 







The influence of mesh size on 

length 2 







the computed mixed mode ratio 

5 mm 

+ 






was negligible for elements 

2 mm 

+ 






lengths between A a = 0.5 mm 

1 mm 

+ 





+ 

(four ply thicknesses) and 








Aa=0.0625 mm (half a ply 

0.5 mm 

+ 





+ 

thickness). Only very long 

0.25 mm 

+ 





+ 

elements (Aa =5 mm ) need to be 

0.125 mm 






+ 

avoided. 

0.0625 mm 






+ 


section 







The mesh refinement study 

length c 3 







showed that only a section of 

1 mm 

+ 






about 1 mm on either side of the 

2 mm 

+ 






delamination front needs to be 

5 mm 

+ 






refined (c=2.0 mm). 

10 mm 

+ 







15 mm 

+ 







20 mm 

+ 







section 







Computed strain energy release 

length d 4 







rates were in good agreement 

5 mm 

+ 5 

+ 

+ 

+ 

+ 

+ 

with the reference solution once 

10 mm 

+ 5 

+ 

+ 

+ 

+ 

+ 

the local zone was extended to 

15 mm 

+ 5 

+ 

+ 

+ 

+ 

+ 

about three times the specimen 








thickness in front and behind the 

20 mm 

+ 

+ 

+ 

+ 

+ 

+ 

delamination front. 

25 mm 

+ 5 

+ 

+ 

+ 

+ 

+ 


30 mm 

+ 5 

+ 

+ 

+ 

+ 

+ 


Contact 


+ 



+ 


MPCs more efficient 

analysis 




+ 


+ 

Contact negligible 


for DCB and ENF specimen c=10 mm, n=20; for SLB specimen c=6 mm, n=24; C3D8I elements used 

2 

the section length, c, was kept constant and the number of elements, n, was varied; c=1 0 mm for DCB specimen, 
c=3 mm for SLB specimen; C3D8I elements used 

3 

the section length, c, was varied as shown; n was modified to keep a constant element length c/n=0.5 mm; C3D8I 
elements used 

4 S4 type shell elements combined with C3D8I solid elements 

5 for DCB-UD24, additionally, a combination of S8R type shell elements with C3D20R solid elements was studied 
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\ / 

Corner 3 Corner 1 

(a) Specimen with crack locations. 


Matrix Crack Branches 


Delamination A 


Initial Matrix Crack 



(b) Comers 1 and 4 


Delamination B1 Delamination B2 



(c) Comers 2 and 3 

Figure 1. Typical damage patterns observed in skin/stringer specimens [21] 
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composite fuselage panel 



shell element model of 
delaminated top laminate 



Figure 2. Application of shell/3D modeling technique to large built-up structures. 
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layup UD24 layup D±30 


B 25.0 mm 25.4 mm 

2h 3.0 mm 4.06 mm 
2L 150.0 mm 150.0 mm 

a 1 1 1 .5 mm 57.2 mm 

P 12.66 N 10.0 N 


(a) Double Cantilever Beam Specimen (DCB) 




layup UD24 

layup 

o 

CO 

+1 

Cl 

B 

25.0 

mm 

25.5 

mm 

2h 

3.0 

mm 

4.06 mm 

2L 

150.0 

mm 

127.0 

mm 

a 

30.0 

mm 

31.8 

mm 

P 

503.0 

N 

100.0 

N 


(b) End Notched Flexure Specimen (ENF) 



layup UD32 

25.0 mm 

2.03 mm 

2.03 mm 
177.8 mm 

34.3 mm 

100.0 N 


layup D±30 

25.4 mm 

2.03 mm 

2.03 mm 
177.8 mm 

34.3 mm 

100.0 N 
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(a) 3D view 



Figure 6. Virtual Crack Closure Technique for comer nodes in twenty noded elements. 




(a) 3D view 



(b) Top view 

Figure 7. Virtual Crack Closure Technique for midside nodes in twenty noded elements. 
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(b) Detail around delamination front (c) Detail around delamination front 

(c= 10 mm, n=20) (c=20 mm, n=40) 

Figure 8. Finite element model of a DCB specimen with UD24 layup 
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(a) Detail around delamination front 
(c=10 mm, n=5) 


(b) Detail around delamination front 
(c= 10 mm, n=20) 



(c) Detail around delamination front 
(c=10 mm, n=40) 

Figure 9. Mesh detail around delamination front 
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Applied Load, N 

Figure 10. Influence of element selection on computed load-displacement 
behavior of a DCB specimen with UD24 layup . 
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Figure 1 1 . Influence of element selection on computed strain energy release rate 
distribution across the width of a DCB specimen with UD24 layup. 
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Figure 12. Influence of refined section on computed strain energy release rate 
distribution across the width of a DCB specimen with UD24 layup. 
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Energy 
Release 
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Figure 13. Influence of number of elements in refined section on computed strain energy 
release rate distribution across the width of a DCB specimen with UD24 layup. 
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(b) Detail (c= d= 5 mm) (c) Detail (d= 10 mm) (d) Detail (d= 15 mm) 






Figure 15. Strain energy release rate distribution across the width of a 
DCB specimen with UD24 layup modeled with 20 noded elements. 


Normalized 
Mode I 
Energy 
Release 
Rate 



Figure 16. Strain energy release rate distribution across the width of a 
DCB specimen with UD24 layup modeled with 8 noded elements. 
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Figure 17. Finite element model of an ENF specimen with UD24 layup 



multi point constraints 
z-direction constrained 

Figure 18. Multi point constraints to prevent contact of delaminated surfaces 


-41 - 



1.500 


i°ooooooo 

jpooooooo 


1.000 


o o o o o o 
X X g Sx>§ Hxg * Xg S >§ s >g 


0.500 


o 

C3D8 

□ 

C3D8I 

o 

C3D8R 

A 

C3D20 

V 

C3D20R 

X 

Kruger [23] 
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Figure 19. Influence of element selection on computed mode II strain energy release rate 
distribution across the width of an ENF specimen with UD24 layup. 
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Figure 20. Influence of element selection on computed mode III strain energy release rate 
distribution across the width of an ENF specimen with UD24 layup. 
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Figure 21. Influence of delamination surface contact on computed strain energy release ratt 
distribution across the width of an ENF specimen with UD24 layup. 
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detail of local 
model around 
delamination front 




Figure 22. Shell/3Dfinite element model of an ENF specimen with UD24 layup 

(c= 5 mm, n= 10, d= 30 mm) 
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Figure 23. Mode II strain energy release rate distribution across the width of an ENF specimen 
with UD24 layup calculated using the shell/3D modeling technique 

0.05 


0.04 


Normalized 0 03 
Mode III 
Energy 
Release 
Rate 0.02 


0.01 


0.00 

-0.4 -0.2 0 0.2 0.4 

y/B 

Figure 24. Mode III strain energy release rate distribution across the width of an ENF specimen 
with UD24 layup calculated using the shell/3D modeling technique. 
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Figure 25. Finite element model of a SLB specimen with UD32 layup. 
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Figure 26. Influence of element selection on computed mode I strain energy release rate 
distribution across the width of a SLB specimen with UD32 layup. 
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Figure 27. Influence of element selection on computed mode II strain energy release rate 
distribution across the width of a SLB specimen with UD32 layup. 
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Figure 28. Influence of element selection on computed mode III strain energy release rate 
distribution across the width of a SLB specimen with UD32 layup. 
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Figure 30. Mode I strain energy release rate distribution across the width of a SLB specimen 
with UD32 layup calculated using the shell/3D modeling technique 



Figure 31. Mode II strain energy release rate distribution across the width of a SLB specimen 
with UD32 layup calculated using the shell/3D modeling technique 
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Figure 32. Mode III strain energy release rate distribution across the width of a SLB specimen 
with UD32 layup calculated using the shell/3D modeling technique 
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Figure 33. Finite element model of SLB specimen with D±30 layup 
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Figure 34. Influence of number of elements in refined section on computed mode I strain energy 
release rate distribution across the width of a SLB specimen with D±3() layup. 
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Figure 35. Influence of number of elements in refined section on computed mode II strain energy 
release rate distribution across the width of a SLB specimen with D±30 layup. 
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Figure 36. Influence of number of elements in refined section on computed mode III strain energy 


release rate distribution across the width of a SLB specimen with D±30 layup. 



Figure 37. Influence of number of elements in refined section on mixed mode ratio 
distribution across the width of a SLB specimen with D±30 layup. 
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Figure 38. Influence of delamination surface contact on computed mode I strain energy 
release rate distribution across the width of SLB specimen with D±30 layup. 
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Figure 39. Influence of delamination surface contact on computed mode II strain energy 
release rate distribution across the width of SLB specimen with D±30 layup. 
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Figure 40. Influence of delamination surface contact on computed mode III strain energy 
release rate distribution across the width of SLB specimen with D±30 layup. 
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Figure 41. Shell/3Dfinite element model of a SLB specimen with D±30 layup 

(c= 3 mm, n=12, d= 20 mm) 
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Figure 42. Mode I strain energy release rate distribution across the width of a SLB specimen 
with D±30 layup calculated using the shell/3 D modeling technique 
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Figure 43. Mode II strain energy release rate distribution across the width of a SLB specimen 
with D±30 layup calculated using the shell/3D modeling technique 
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Figure 44. Mode III strain energy release rate distribution across the width of a SLB specimen 
with D±30 layup calculated using the shell/3D modeling technique 
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Figure 45. Finite element model of a DCB specimen with D±30 layup 
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Figure 46. Influence of delamination surface contact on computed strain energy release rate 
distribution across the width of a DCB specimen with D±30 layup modeled with 8 noded elements 
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Figure 47. Shell/3Dfinite element model of a DCB specimen with D±30 layup 
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Figure 48. Mode I strain energy release rate distribution across the width of a 
DCB specimen with D±30 layup calculated using the shell/3D modeling technique 
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Figure 50. Influence of delamination surface contact on computed mode II strain energy 
release rate distribution across the width of an ENF specimen with D±30 layup. 
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Figure 51. Influence of delamination surface contact on computed mode III strain energy 
release rate distribution across the width of an ENF specimen with D±30 layup. 
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Figure 52. Shell/3Dfinite element model of an ENF specimen with D±30 layup 

(c= 3 mm, n=12, d= 30 mm) 
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Figure 53. Mode II strain energy release rate distribution across the width of an ENF specimen 
with D±30 layup calculated using the shell/3D modeling technique 
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Figure 54. Mode III strain energy release rate distribution across the width of an ENF specimen 
with D±30 layup calculated using the shell/3D modeling technique 
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